## implementation of Next Subvolume Method [6]

source('subvolume_queue.R')

spatial.setup <- function(ctx, X0, transitions, V, W, H, D)
{
	ctx$transitions <- transitions
	ctx$V <- V
	ctx$W <- W
	ctx$H <- H
	ctx$D <- D
	ctx$X <- X0  # step 2

	with(ctx, {
		# step 1

		Nmatrix <- matrix(rep(1:(W*H),4), ncol=4)   # connectivity/adjacency matrix
		Nmatrix[(1:(W*H))[-seq(1,W*H,W)],1] <- Nmatrix[(1:(W*H))[-seq(1,W*H,W)],1]-1
		Nmatrix[(1:(W*H))[-seq(W,W*H,W)],2] <- Nmatrix[(1:(W*H))[-seq(W,W*H,W)],2]+1
		if(W > 1)
			Nmatrix[(W+1):(W*H),3] <- Nmatrix[(W+1):(W*H),3]-W
		if(H > 1)
			Nmatrix[1:(W*(H-1)),4] <- Nmatrix[1:(W*(H-1)),4]+W
		N <- apply(Nmatrix != matrix(rep(1:(W*H),4),ncol=4), 1, sum)

		# step 3,4,5

		R <- apply(X, 1, function(x) sum(transitions(x)))
		S <- N*c(X %*% D)
		tao <- -log(runif(W*H))/(R+S)

		# step 6

		Q <- matrix(c(tao, 1:(W*H)), nrow=W*H)  # queue array "binary tree"
		Qarray <- 1:(W*H)                       # reverse lookup table
		lapply(1:(W*H), function(x) sort.queue(ctx, Qarray[x]))  # ugly: could be faster

		t <- 0
	})
}

spatial.step <- function(ctx)
{
	with(ctx, {
		# step 7

		lambda <- Q[1,2]
		t <- tao[lambda]
		rand <- runif(1)
		RS <- R[lambda]+S[lambda]
		if(RS == 0)  # bug in the algorithm: if no molecules, div/0
			return()
		reaction <- R[lambda]/RS

		if(rand < reaction) {
			# step 8: reaction event
			# a)
			rand <- rand/reaction
			# b)
			r <- cumsum(transitions(X[lambda,]))
			i <- which(rand < r/R[lambda])[1]  # index reaction
			X[lambda,] <- X[lambda,] + V[i,]
			# c)
			rand <- runif(1)
			R[lambda] <- sum(transitions(X[lambda,]))
			S[lambda] <- N[lambda]*sum(X[lambda,]*D)
			tao[lambda] <- -log(runif(1))/(R[lambda]+S[lambda]) + t
			# d)
			Q[1,1] <- tao[lambda]
			sort.queue(ctx, 1)
		}
		else {
			# step 9: diffusion event
			# a)
			rand <- (rand-reaction)/(1-reaction)
			d <- D*X[lambda,]
			i <- which(runif(1) < cumsum(d)/sum(d))[1]  # index molecule
			# b)
			dir <- trunc(rand*4)+1
			gamma <- Nmatrix[lambda,dir]
			# c)
			X[lambda,i] <- X[lambda,i] - 1
			X[gamma,i] <- X[gamma,i] + 1
			# d)
			R[lambda] <- sum(transitions(X[lambda,]))
			S[lambda] <- N[lambda]*sum(X[lambda,]*D)
			tao[lambda] <- -log(runif(1))/(R[lambda]+S[lambda]) + t
			R[gamma] <- sum(transitions(X[gamma,]))
			S[gamma] <- N[gamma]*sum(X[gamma,]*D)
			tao[gamma] <- -log(runif(1))/(R[gamma]+S[gamma]) + t
			# e)
			Q[1,1] <- tao[lambda]
			sort.queue(ctx, 1)
			Q[Qarray[gamma],1] <- tao[gamma]
			sort.queue(ctx, Qarray[gamma])
		}
	})
}

## References:
# 1. lotka-volterra model parameters from example at
#    https://darrenjw.wordpress.com/tag/gillespie/
# 2. R image documentation
#    http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/image.html
# 3. nice plot
#    http://stackoverflow.com/questions/12918367/in-r-how-to-plot-with-a-png-as-background
# 4. Higham, D. J. (2008). Modeling and simulating chemical reactions. SIAM
#    review, 50(2), 347-368.
# 5. Banks, H. T., Hu, S., Joyner, M., Broido, A., Canter, B., Gayvert, K., &
#    Link, K. (2012). A comparison of computational efficiencies of stochastic
#    algorithms in terms of two infection models. Mathematical biosciences and
#    engineering: MBE, 9(3), 487-526.
# 6. Elf and Ehrenberg (2004) Spontaneous separation of bi-stable biochemical
#    systems into spatial domains of opposite phases. Systems Biolog.
#    Supplementary material: Next Subvolume Method algorithm
